1. DADES DE MALLORCA
Maria Àngels Jaume
2023-02-04
1. LECTURA DE DADES
Un cop descarregades les dades, les llegim i les grafiquem en forma de sèrie temporal. Visualitzam la sèrie de Mallorca, Menorca, Eivissa i Formentera, i les Illes Balears en una mateixa gràfica:
turisme <- read_excel("IBESTAT.xls")
turisme$Data <- gsub("M","-",turisme$Data)
gastos.ts<-ts(turisme[-1], start=c(2015,10), frequency = 12)plot(gastos.ts, plot.type = "single", col = 1:ncol(gastos.ts), ylim=c(0,1300), xlab="Any", ylab="Despeses mensuals en €", main='Series de cada illa')
legend("bottomright", colnames(gastos.ts), col=1:ncol(gastos.ts), lty=1.5, cex=.85)
Les dades van des d’octubre de 2015 a novembre de 2022 (llavors tenim 6 cicles complets). Podem observar que totes les illes segueixen més o menys un mateix patró (amb una diferència entre les series de Menorca i Eivissa i Formentera davant les de Mallorca i les Illes Balears) i l’efecte de la COVID és veu clarament en 2020. També observam una clara estacionalitat que pareix tenir els seus valors més alts en estiu i els més baixos en hivern.
En aquest document estudiarem l’illa de Mallorca, diguem
mall al conjunt de dades que fan referència a l’illa de
Mallorca.
mall <- data.frame(x=1:86, y=turisme$Mall)
head(mall)## x y
## 1 1 895.45
## 2 2 850.73
## 3 3 834.39
## 4 4 773.57
## 5 5 858.62
## 6 6 814.16La variable \(X\) ens indica el mes en que ens trobam (hem de tenir en compte que les dades comencen el 10 de novembre 2015, llavors el valor 1 correspon a aquest més i segueix fins al mes 86, que correspon a novembre de 2022), mentre que la variable \(Y\) ens indica les despeses dels turistes en aquells mes en euros.
Abans d’aplicar algun model i posar-nos a analitzar les dades, visualitzem-les per tenir una idea de qué ens podem trobar a l’hora d’aplicar els models:
serie_mall <- ts(mall$y,start=c(2015,10),frequency = 12)
plot.ts(serie_mall, main="Illa de Mallorca", ylab="Despeses mensuals en €", xlab="Any")
Per veure millor l’estacionalitat de la sèrie visualitzem cada període mensualment:
seasonplot(serie_mall, col=c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),year.labels=TRUE, main="Estacionalitat de Mallorca", xlab="Mes", ylab=" Despeses en €")
legend("bottomright",
legend = c(2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
fill = c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),
border = "black") 
Observam que efectivament, així com es comenta en el treball, el període de confinament va ser de principi de març, on comença a decréixer el valor dels preus, fins a finals de juny, on pareix que s’han tornat a recuperar els valors anteriors a la COVID. A més, també observam que les temporades d’estiu són més altes que les de l’hivern.
Dividirem la sèrie en tres trams:
Des de 10-2015 fins 12-2018 (serie1_mall), amb el que predirem un període. És a dir, predirem de 1-2019 a 12-2019. (amb aquest comprovam que els mètodes funcionen)
Des de 10-2015 fins 12-2019 (serie2_mall), amb la que predirem el període de la COVID (és a dir, l’any 2020)
Des de 10-2015 fins 12-2020 (serie3_mall), amb la que predirem el període després de la COVID (és a dir, l’any 2021)
Visualitzem-los:
serie1_mall <- ts(mall$y[1:39],start=c(2015,10),frequency = 12)
serie2_mall <- ts(mall$y[1:51],start=c(2015,10),frequency = 12)
serie3_mall <- ts(mall$y[1:63],start=c(2015,10),frequency = 12)
par(mfrow=c(2,2), mar=c(4,4,4,1)+.1)
plot.ts(serie1_mall, main="Sèrie abans de la COVID \nmenys un cicle" , xlab="Any", ylab="despeses mensuals")
plot.ts(serie2_mall, main="Sèrie abans de la COVID", xlab="Any", ylab="despeses mensuals")
plot.ts(serie3_mall, main="Sèrie abans i durant de la COVID", xlab="Any", ylab="despeses mensuals")
plot.ts(serie_mall, main="Sèrie completa", xlab="Any", ylab="despeses mensuals")
2. AJUST DELS MODELS:
Abans d’aplicar algun model, estudiem un poc el tram a analitzar:
serie_mall.lm <- lm(y~x, data=mall)
plot.ts(mall$y, main = "Mallorca ", ylab = "Despeses menusals en €", xlab="Índex de cada mes")
#dibuixam la recta de regressió sobre els punts
abline(serie_mall.lm, col='red')
Si dibuixam la recta de regressió sobre les nostres dades, tot i que aquestes estan molt disperses i no s’ajusten bé a la recta, podem observar que la tendència és una mica creixent, encara que es manté més o menys constant.
boxplot(serie_mall~cycle(serie_mall), xlab = "Mes", ylab = "Despeses en €", main="Boxplot de Mallorca")
Podem observar també la presència d’estacionalitat, que prenen els seus valors màxims a la temporada d’estiu i els seus mínims en l’hivern, fet que correspon amb les temporades turístiques a l’illa.
Aplicarem diversos models per ajustar la nostra sèrie i fer-ne una predicció per llavors comparar quin és el millor.
Vegem ara com actua cada model:
2.1. ABANS DE LA COVID
- REGRESSIÓ SEGMENTADA
Com hem comentat anteriorment la recta de regressió no s’ajusta bé a les dades, de fet el valor del \(R^2\) és molt baix:
serie1_mall.lm <- lm(y~x, data=mall[1:39,])
summary(serie1_mall.lm)$adj.r.squared## [1] 0.02188284Per això utilitzam el paquet segmented per ajustar a una
recta de regressió a trossos.
Anem a utilitzar la funció selgmented() per veure quants
de punts de canvi detecta:
punts_canvi_serie1_mall <-selgmented(serie1_mall.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 ..
##
## BIC to detect no. of breakpoints:
## 0 1 2 3 4 5 6 7
## 370.5453 374.5279 381.6098 376.5556 372.6382 357.6041 351.5419 337.9820
## 8 9 10
## 344.0487 349.0722 355.4233
##
## No. of selected breakpoints: 7Aplicam la funció segmented() que ens calcula la
regressió segmentada:
serie1_mall.seg <- segmented(serie1_mall.lm, seg.Z=~x, psi = c(3,10,17,23,28,35))Aquesta funció ens demana que introduïm els valors on es troben els punts de canvi, i aquesta ens dóna el valor estimat. Aquests punts de canvi són:
summary(serie1_mall.seg)$psi## Initial Est. St.Err
## psi1.x 3 6.557347 0.7808007
## psi2.x 10 10.579679 0.4882573
## psi3.x 17 16.351485 0.4500235
## psi4.x 23 22.934881 0.4249014
## psi5.x 28 27.781744 0.4721648
## psi6.x 35 34.980655 0.4401346Que corresponen, aproximadament, a: 4-2016, 8-2016, 1-2017, 8-2017, 1-2018, 8-2018. Llavors tenim que els punts més alts es donen per agost i els més baixos per gener.
Obtenim que el valor de \(R^2\) per a la regressió segmentada és bastant alt
summary(serie1_mall.seg)$adj.r.squared## [1] 0.8063167Anem a visualitzar la regressió segmentada damunt les nostres dades
#Per graficar-ho
p_serie1_mall <- ggplot(mall[1:39,], aes(x=x, y=y), color="black") +
geom_line()+ #dibuixam les nostres dades en línia contínua
ggtitle('Regressió lineal i segmentada sobre les dades de Mallorca \nabans de la COVID') +
xlab('Índex del mes') +
ylab('Despeses mensuals en €') +
ylim(c(0,1300))
my.coef1_mall <- coef(serie1_mall.lm) #coeficients de la recta de regressió lineal
p_serie1_mall <- p_serie1_mall + geom_abline(intercept=my.coef1_mall[1], slope=my.coef1_mall[2], color="green") #afegim la recta de regressió lineal
my.model1_mall <- data.frame(x=1:39, y=fitted(serie1_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie1_mall <- p_serie1_mall + geom_line(data=my.model1_mall, aes(x=x,y=y), color="red") #afegim la recta de regressió segmentada
my.lines1_mall <- serie1_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie1_mall <- p_serie1_mall + geom_vline(xintercept=my.lines1_mall, linetype="dashed") #línies verticals en els punts de canvi
p_serie1_mall <- p_serie1_mall + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
ggplotly(p_serie1_mall)Visualment també es pot observar que la recta de regressió segmentada s’ajusta millor a les nostres dades.
Anem a calcular les equacions d’aquestes rectes, sabem que les rectes tenen la forma \(y=mx+n\), on \(m\) és la pendent i \(n\) el valor de tall de l’eix y.
Hi ha una funció del paquet segmented que ens dóna
aquests valors de \(n\):
intercept(serie1_mall.seg)## $x
## Est.
## intercept1 882.18
## intercept2 354.63
## intercept3 1740.80
## intercept4 -300.59
## intercept5 2844.30
## intercept6 -606.28
## intercept7 3780.00També tenim una altra funció que ens calcula les pendents:
slope(serie1_mall.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -12.674 11.4980 -1.1023 -36.355 11.006
## slope2 67.777 21.5110 3.1509 23.475 112.080
## slope3 -63.248 11.4980 -5.5008 -86.928 -39.567
## slope4 61.599 11.4980 5.3574 37.919 85.279
## slope5 -75.525 15.2100 -4.9654 -106.850 -44.199
## slope6 48.679 9.0899 5.3553 29.958 67.400
## slope7 -76.713 15.2100 -5.0435 -108.040 -45.387Aleshores, la nostra recta segmentada és:
\(\hat{y}= \left\{ \begin{array}{lcc} -12.674x + 882.18 & si & x \leq \psi_1 \\ \\ 67.777x + 354.63 & si & \psi_1 < x \leq \psi_2 \\ \\ -63.248x + 1740.8 & si & \psi_2 < x \leq \psi_3 \\ \\ 61.599x - 300.59 & si & \psi_3 < x \leq \psi_4 \\ \\ -75.525x + 2844.3 & si & \psi_4 < x \leq \psi_5 \\ \\ 48.679x - 606.28 & si & \psi_5 < x \leq \psi_6 \\ \\ -76.713x + 3780 & si & \psi_6 < x \\ \end{array} \right.\)
on \(\psi_1= 6.56\), \(\psi_2= 10.58\), \(\psi_3= 16.35\), \(\psi_4= 22.93\), \(\psi_5= 27.78\) i \(\psi_6= 34.98\)
Vegem els errors del model (ens interessa el RMSE):
accuracy(serie1_mall.seg)## ME RMSE MAE MPE MAPE MASE
## Training set 2.915581e-15 38.51009 32.1737 -0.1792044 3.633043 0.3559833Anem a fer la predicció de 1-2019 a 12-2019.
El darrer punt de canvi que tenim és en agost de 2018 i sabem, seguint el patró obtingut amb les dades d’entrenament, que els següents és donaran en gener de 2019, agost de 2019 i gener de 2020. Necessitam calcular les pendents de les rectes d’entre agost de 2018 i gener de 2020, per poder fer la predicció i calcular l’error.
Per calcular les pendents de les noves rectes farem la mitjana de les pendents anteriors. De les pendents ja calculades en el model obviarem la primera i la darrera, ja que no són vàlides. Per tant,
- La pendent de 8-2018 a 1-2019 i de 8-2019 a 1-2020 serà : (-63.248 - 75.525)/2 = -69.387
- La pendent de 1-2019 a 8-2019 serà : (66.777 + 61.599 + 48.679)/3 =59.018
Llavors els nous segments predits seran
\(\hat{y}= \left\{ \begin{array}{lcc} -69.387x + n_1 & si & x \leq \psi_7 \\ \\ 59.018x + n_2 & si & \psi_7 < x \leq \psi_8 \\ \\ -69.387x + n_3 & si & \psi_8 < x \\ \end{array} \right.\)
on \(\psi_7 = 40\) i \(\psi_8 = 47\).
Calculem aquests \(n_i\), \(i\in \{1,2,3\}\):
La darrera recta vàlida de la regressió segmentada és \(y = 48.679x - 606.28\)
Per \(n_1\);
Si \(x = 34.98\), que és el darrer punt de canvi, la \(y = 48.679 \cdot 34.98 - 606.28 = 1096.511\).Per tant, el punt d’intersecció de la primera recta a predir és \(1096.511 = -69.387 \cdot 34.98 + n_1\). Llavors \(n_1 = 3523.67\).
Per \(n_2\);
Si \(x= 40\), que és el primer punt de canvi predit, tenim \(y = -69.387 \cdot 40 + 3523.67 = 748.19\). Per tant, el segon punt d’intersecció amb l’eix OY és \(748.19 = 59.018 \cdot 40 + n_2\), \(n_2 = -1612.53\).
Per \(n_3\);
Si \(x=47\), que és el segon punt de canvi predit, \(y = 59.018 \cdot 47 -1612.53 = 1161.315\). Per tant, el tercer punt d’intersecció amb l’eix OY és \(1161.315 = -69.387 \cdot 47 + n_3\), \(n_3 = 4422.48\).
Aleshores, la predicció per al 2019 amb la regressió segmentada és
\(\hat{y}= \left\{ \begin{array}{lcc} -69.387x + 3523.67 & si & x \leq \psi_7 \\ \\ 59.018x -1612.53 & si & \psi_7 < x \leq \psi_8 \\ \\ -69.387x + 4422.48 & si & \psi_8 < x \\ \end{array} \right.\)
on \(\psi_7 = 40\) i \(\psi_8 = 47\).
#Per graficar-ho
p1_serie_mall <- ggplot(mall[1:51,], aes(x=x, y=y)) + geom_line()+
ylim(c(500,1300)) +#dibuixam les nostres dades en línia contínua
ggtitle('Predicció de Mallorca amb el model \nde regressió segmentada abans de la COVID') +
xlab('Índex del mes') +
ylab('Despeses mensuals en €')
my.model1_mall <- data.frame(x=1:39, y=fitted(serie1_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada
p1_serie_mall <- p1_serie_mall + geom_line(data=my.model1_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines1_mall <- serie1_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats
p1_serie_mall <- p1_serie_mall + geom_vline(xintercept=my.lines1_mall, linetype="dashed") #línies verticals en els punts de canvi
p1_serie_mall <- p1_serie_mall+ geom_vline(xintercept=0) + geom_hline(yintercept=500) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
p1_serie_mall <- p1_serie_mall + geom_abline(intercept = 3523.67, slope=-69.387, colour="green") +
geom_abline(intercept = -1612.53, slope=59.018, colour="blue") +
geom_abline(intercept = 4422.48, slope=-69.387, colour="orange")
ggplotly(p1_serie_mall)Observam que les noves rectes predites s’ajusten més o menys bé a la sèrie.
Calculem l’error de la predicció:
o1_mall<-c(serie_mall[40:51]) #observacions reals de la sèrie de gener de 2019 a desembre de 2019
v1=c(40:47)
f1 <- sapply(v1, function(x) 59.018*x-1612.53) #predicció de gener 2019 a agost 2019
v2=c(48:51)
f2 <- sapply(v2, function(x) -69.387*x + 4422.48)#predicció de setembre 2019 a desembre 2019
p1_mall<- c(f1,f2) #predicció de gener de 2019 a desembre de 2019
sqrt(sum((p1_mall-o1_mall)^2)/12) #error de la predicció## [1] 52.35583Aleshores aquest model pareix ser bastant bo, ja que per fer la predicció un any endavant només tenim un error d’uns 43 euros.
- DECOMPOSE
Recordem la sèrie
plot.ts(serie1_mall, main= "Mallorca abans de la COVID", xlab="Any", ylab="Despeses mensuals en €")
Ja hem dit anteriorment que té una tendència mínima i podem observar, també, que no hi ha variabilitat. Llavors podem suposar un model additiu.
El mètode clàssic de descomposició separa la sèrie en subseries corresponents a la tendència, la estacionalitat i el renou.
Aquestes components es poden combinar de manera additiva o multiplicativa. En el nostre cas utilitzam el model additiu: \(y_t = \mu_t + S_t + a_t\)
decompose_s1_mall<-decompose(serie1_mall)
plot(decompose_s1_mall, xlab="Any")
El decompose, per calcular aquestes noves tendències, el
que fa és agafar les sis tendències anteriors i les sis posteriors de la
sèrie original i en fa la mitjana. És per això que la primera que
obtenim és en abril de 2016 i la darrera en juny de 2018. Notem que per
a la predicció ens quedarem amb el valor de la darrera tendència del
decompose.
t1_mall <- decompose_s1_mall$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t1_mall## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 899.8458 899.0708 896.9637 895.5496 888.9975
## 2017 891.8125 896.2783 901.1758 905.1896 907.7217 908.4975 909.7354 914.6108
## 2018 924.3500 924.6162 923.1429 921.8554 921.3479 920.0171 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 885.1633 886.3112 887.5750 888.9183
## 2017 920.9775 923.1212 922.9363 923.5758
## 2018 NA NA NA NAEls valors de les components estacionals els calcula fent la mitjana per mesos, és a dir, per calcular la componen de gener, agafa els valors de tots els geners que tenim a la sèrie original i en fa la mitjana. Per tant, només tenim 12 valors, un per a cada mes.
s1_mall <- decompose_s1_mall$seasonal #components estacionals
s1_mall <- s1_mall[4:15] #estacionalitat de gener a desembre
s1_mall## [1] -174.312598 -107.883640 -54.675723 -66.333293 6.116846 47.367541
## [7] 169.981152 195.059485 82.053235 67.382402 -116.501973 -48.253432Anem a fer la predicció d’aquest model
a1_mall <- c(s1_mall[7:12],s1_mall) #estacionalitat de juliol 2018 a desembre 2019 (és per poder fer la predicció)
pred1_decompose_mall <- sapply(a1_mall, function(x) 920.0171 + x) #predicció de juliol de 2018 a desembre 2019 (el valor 920.0171 és la darrera tendència obtinguda amb el decompose)p1_dec_mall<-c(serie1_mall[1:33], pred1_decompose_mall) #valors de la sèrie original més els predits
A1_mall<- data.frame("x" = mall[1:51,]$x, "y" = mall[1:51,]$y, "p"= p1_dec_mall) #cream un data frame on hi trobam una columna amb els mesos, x, una amb els valors originals de la sèrie, y, i una altra amb les prediccions, p.
#ho dibuxam
grafica1_mall_dec <- ggplot(A1_mall)+
ylim(c(600,1200)) +
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+ geom_vline(xintercept=0) +
geom_hline(yintercept=600)+
labs(title="Predicció de Mallorca abans de la COVID amb el model de descomposició", x="Índex del mes", y="Despeses mensuals en €")+
theme(panel.background = element_blank())
grafica1_mall_dec
Podem observar que aquest model fa una predicció bastant bona. No obstant això observam que, com el quart cicle és un poc més alt que els anteriors, la predicció no arriba a captar els valors més alts del darrer cicle.
Calculem l’error de la predicció:
(Notem que la predicció és des de juliol de 2018 a desembre de 2019 i per calcular l’error només volem de gener de 2019 a desembre de 2019.)
sqrt(sum((c(serie_mall[40:51]- pred1_decompose_mall[7:18]))^2)/12) #error predicció de gener a desembre de 2019## [1] 32.40833Aleshores l’error que és comet amb el model clàssic de descomposició és de 32 euros, aproximadament.
Notem que també tenim una altra instrucció a R per fer prediccions
d’una sèrie, la funció predict(). Aquesta és basa en un
model ETS, anem a veure la predicció:
prediccio1_mall <- predict(serie1_mall,h=12)
plot(prediccio1_mall, xlab="Any", ylab="Despeses mensuals en €")
Pareix que la predicció és bastant bona, ja que el cicle predit segueix un mateix patró que els anteriors. Anem a veure la predicció juntament amb la sèrie original:
df_prediccio1_mall <- data.frame(prediccio1_mall)
p1_ets_mall <- data.frame("x"= 1:51, "PointForecast"=c(serie1_mall,df_prediccio1_mall$Point.Forecast), "Lo80" = c(rep(NA,39),df_prediccio1_mall$Lo.80), "Hi80" = c(rep(NA,39),df_prediccio1_mall$Hi.80), "Lo95" = c(rep(NA,39),df_prediccio1_mall$Lo.95),"Hi95" = c(rep(NA,39),df_prediccio1_mall$Hi.95))
grafica_pred1_ets <- ggplot((mall[1:51,]))+ ylim(c(600,1200))+
geom_ribbon(data = p1_ets_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p1_ets_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p1_ets_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(A,N,A) a Mallorca abans de ls COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=600)+ theme(panel.background = element_blank())
grafica_pred1_ets
Podem observar que la predicció és bastant bona, ja que continua seguint un mateix patró. I l’error comés seria de 37 euros.
accuracy(prediccio1_mall,serie2_mall)## ME RMSE MAE MPE MAPE MASE
## Training set 1.886318 26.85308 22.27401 0.1691627 2.499489 0.6617645
## Test set 19.065194 37.02366 31.47709 1.6792842 3.289203 0.9351895
## ACF1 Theil's U
## Training set -0.01565369 NA
## Test set 0.43720293 0.4375391- SARIMA
Anem a veure quin model obtenim considerant un model estacional.
Pel que hem vist anteriorment, podem considerar que no hi ha tendència, llavors no fa falta fer cap diferència a la part regular, no obstant, sí que cal fer una diferenciació d’orde 12 per l’estacionalitat. Vegem l’ACF i el PACF de la sèrie a predir:
par(mfrow=c(1,2))
acf(serie1_mall)
pacf(serie1_mall)
L’ACF ens està indicant un model de mesures mòbils d’ordre 2, MA(2). És a dir, té en compte els errors comesos en els dos mesos anteriors a l’hora de fer la predicció. A més, el PACF ens indica un model autoregressiu d’ordre 1, AR(1). Aquest fet ens està indicant que a l’hora de fer la predicció es té en compte el valor obtingut al mes anterior.
Llavors, per a la part regular obtenim un ARIMA(1,0,2).
Fem una diferenciació a la part estacional, és a dir, d’ordre 12
serie1_mall_dif <- diff(serie1_mall,12)
plot(serie1_mall_dif, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")
Aquesta és la sèrie sense estacionalitat, vegem com es modifiquen l’ACF i el PACF.
par(mfrow=c(2,1))
acf(serie1_mall_dif, lag=36)
pacf(serie1_mall_dif,lag=36)
Observam que no hi ha barres significatives, llavors no tenim cap model autoregressiu ni de mesures mòbils a la part estacional. És a dir, no es tenen en compte els valors obtinguts dels cicles anteriors.
Per tant hem obtingut un ARIMA(1,0,2)(0,1,0)[12].
Vegem quin model els proposa R:
auto.arima(serie1_mall)## Series: serie1_mall
## ARIMA(0,0,0)(1,1,0)[12] with drift
##
## Coefficients:
## sar1 drift
## -0.5576 0.8407
## s.e. 0.1924 0.4397
##
## sigma^2 = 1342: log likelihood = -136.72
## AIC=279.45 AICc=280.49 BIC=283.33R ens proposa un ARIMA(0,0,0)(1,1,0)[12], és a dir només fa una diferenciació a la part estacional i considera que per predir un mes s’han de tenir en compte els errors d’aquest mateix mes però en el cicle anterior.
Per tant les propostes de model ARIMA són:
model1_mall<-arima(serie1_mall, order=c(1,0,2), seasonal = list(order=c(0,1,0), period=12))
model2_mall <- arima(serie1_mall, order=c(0,0,0), seasonal = list( order=c(1,1,0), period=12))Els errors d’aquests dos models obtinguts són
accuracy(model1_mall)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 4.916634 32.66904 22.17796 0.380845 2.549795 0.3062364 -0.03442627accuracy(model2_mall)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 8.180359 32.10763 21.18521 0.6723209 2.372129 0.2925283 0.2063501Observam que no hi ha molta diferència d’un model a un altre.
Anem a fer la predicció de 2019 d’aquests models :
#ens torna un data frame amb els punts predits, un interval del 80% i un altre del 95%
fc_m1_mall <-forecast(model1_mall, h=12)
fc_m2_mall <- forecast(model2_mall,h=12)Visualitzem-les:
#MODEL 1
pre1 <- data.frame("Point Forecast" = serie1_mall, "Lo 80" = rep(NA,39), "Hi 80"= rep(NA,39), "Lo 95" = rep(NA,39), "Hi 95" = rep(NA,39)) #dades abans de la predicció amb NA als intervals ja que només els volem per la predicció
pred_m1_mall <-data.frame(fc_m1_mall)
#vectors amb les dades de la sèrie i la predicció amb els intervals
a1_mall <- c(pre1$Point.Forecast,pred_m1_mall$Point.Forecast)
b1_mall <- c(pre1$Lo.80, pred_m1_mall$Lo.80)
c1_mall <- c(pre1$Hi.80, pred_m1_mall$Hi.80)
d1_mall <- c(pre1$Lo.95, pred_m1_mall$Lo.95)
e1_mall <- c(pre1$Hi.95, pred_m1_mall$Hi.95)
grafica_m1_mall <- data.frame("x" = 1:51, "PointForecast" = a1_mall, "Lo80" = b1_mall, "Hi80"= c1_mall, "Lo95" = d1_mall, "Hi95" = e1_mall)
grafica1_mall <- ggplot(mall[1:51,])+ ylim(c(600,1200))+
geom_ribbon(data = grafica_m1_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m1_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m1_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció de Mallorca amb el model ARIMA (1,0,2)(0,1,0)[12] abans de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=600)+ theme(panel.background = element_blank())
grafica1_mall
Podem observar que la predicció (línia blava) és bastant bona ja que és molt propera a la de les dades reals (la vermella). A més, les bandes dels intervals de confiança no són molt amplis, el que ens indica que la probabilitat d’error no és molt gran.
#MODEL 2
pred_m2_mall <-data.frame(fc_m2_mall)
a2_mall <- c(pre1$Point.Forecast,pred_m2_mall$Point.Forecast)
b2_mall <- c(pre1$Lo.80, pred_m2_mall$Lo.80)
c2_mall <- c(pre1$Hi.80, pred_m2_mall$Hi.80)
d2_mall <- c(pre1$Lo.95, pred_m2_mall$Lo.95)
e2_mall <- c(pre1$Hi.95, pred_m2_mall$Hi.95)
grafica_m2_mall <- data.frame("x" = 1:51, "PointForecast" = a2_mall, "Lo80" = b2_mall, "Hi80"= c2_mall, "Lo95" = d2_mall, "Hi95" = e2_mall)
grafica2_mall <- ggplot(mall[1:51,])+ ylim(c(600,1200))+
geom_ribbon(data = grafica_m2_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m2_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m2_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció de Mallorca amb el model ARIMA (0,0,0)(1,1,0)[12] abans de la COVID", x="índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=600)+ theme(panel.background = element_blank())
grafica2_mall
Amb aquest model la predicció pareix un poc més bona, ja que les bandes blaves són més fines i la línia blava, que són els valors predits, s’ajusten millor als reals (els vermells).
Vegem quins són els errors de la predicció:
accuracy(fc_m1_mall, serie_mall[40:51], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set 4.916634 32.66904 22.17796 0.380845 2.549795 0.3062364
## Test set 13.314448 39.50200 32.06322 1.095781 3.483925 0.4427334
## ACF1
## Training set -0.03442627
## Test set NAaccuracy(fc_m2_mall,serie_mall[40:51], h=12)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 8.180359 32.10763 21.18521 0.6723209 2.372129 0.2925283 0.2063501
## Test set 17.204294 30.21333 24.52273 1.6109194 2.569169 0.3386133 NAEfectivament, com hem vist que pareixia amb les gràfiques, la predicció del segon model dona menys error (30.21) que el primer(39.5), però la diferència és mínima.
Els errors o residus són útils per comprovar si un model ha capturat adequadament la informació de les dades. Un bon mètode de pronòstic produirà residus amb les següents propietats:
Els residus no estan correlacionats. Si hi ha correlació, queda estructura per eliminar.
Els residus tenen mitjana zero.
La variància és constant.
Es distribueixen normalment.
Estudiem-los:
checkresiduals(fc_m1_mall)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2)(0,1,0)[12]
## Q* = 3.4598, df = 5, p-value = 0.6295
##
## Model df: 3. Total lags used: 8checkresiduals(fc_m2_mall)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,0)(1,1,0)[12]
## Q* = 6.5782, df = 7, p-value = 0.4741
##
## Model df: 1. Total lags used: 8Amb les gràfiques del checkresiduals podem observar que,
en els dos models, la sèrie dels errors no mostra cap estructura i les
barres dels ACF es troben dins les bandes de confiança. A més, s’observa
que les distribucions dels errors no pareixen seguir una normal. Aquesta
funció també ens proporciona el p-valor del test de Ljung-Box, que
avalua si els residus del model s’ajusten a un procés de renou blanc, és
a dir, si els residus són independents i no correlacionats. El model
ARIMA(1,0,2)(0,1,0)[12] pareix seguir millor l’estructura de la
gaussiana.
Vegem els qqPlot i apliquem tests de normalitat:
par(mfrow=c(1,2))
qqPlot(fc_m1_mall$residuals, id=FALSE, mean=mean(fc_m1_mall$residuals), sd=sd(fc_m1_mall$residuals))
qqPlot(fc_m2_mall$residuals, id=FALSE, mean=mean(fc_m2_mall$residuals), sd=sd(fc_m2_mall$residuals))
Podem observar que en el primer model tenim més valors centrats al 0 i menys valors fora de la banda. Tot i així els valors no es troben sobre la diagonal, fet que ens indica que no hi ha normalitat.
shapiro.test(fc_m1_mall$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m1_mall$residuals
## W = 0.91601, p-value = 0.006562shapiro.test(fc_m2_mall$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m2_mall$residuals
## W = 0.89116, p-value = 0.001243lillie.test(fc_m1_mall$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m1_mall$residuals
## D = 0.17804, p-value = 0.003142lillie.test(fc_m2_mall$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m2_mall$residuals
## D = 0.25238, p-value = 1.126e-06Tots els tests ens donen un valor molt baix, concretament menor a 0,05. Fet que ens indica que no hem d’acceptar la hipòtesi de que la distribució segueix una normal.
Vegem un resum dels errors que cometen cada un dels models
anteriors:
| reg. segmentada | descomposició clàssica | ETS(A,N,A) | ARIMA (1,0,2)(0,1,0)[12] | ARIMA (0,0,0)(1,1,0)[12] | |
|---|---|---|---|---|---|
| error model | 38.51 | 26.85 | 32.6690 | 32.1076 | |
| error predicció | 52.356 | 32.40833 | 37.02366 | 39.5020 | 30.21333 |
2.2. DURANT LA COVID
- REGRESSIÓ SEGMENTADA
Ara el valor de \(R^2\) per a la regressió lineal és:
serie2_mall.lm <- lm(y~x, data= mall[1:51,])
summary(serie2_mall.lm)$adj.r.squared## [1] 0.02816003És a dir, un valor molt baix. Llavors aplicam la regressió segmentada seguint el mateix procediment que abans.
Vegem quants de punts de canvi tenim i identifiquem-los:
punts_canvi_serie2_mall <-selgmented(serie2_mall.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 ..
##
## BIC to detect no. of breakpoints:
## 0 1 2 <NA> 3 <NA> 4 5
## 489.3701 496.5603 503.9916 504.9916 505.6789 506.6789 500.9452 502.8412
## 6 7 8
## 496.1681 453.0991 440.6423
##
## No. of selected breakpoints: 8serie2_mall.seg <- segmented(serie2_mall.lm, seg.Z=~x, psi = c(4,10,17,23,28,35,40,47))
summary(serie2_mall.seg)$psi## Initial Est. St.Err
## psi1.x 4 4.263718 0.7422268
## psi2.x 10 11.000003 0.5850117
## psi3.x 17 16.430363 0.4682119
## psi4.x 23 22.901107 0.4187605
## psi5.x 28 27.805715 0.4687236
## psi6.x 35 34.907009 0.4340051
## psi7.x 40 40.227949 0.3993340
## psi8.x 47 46.922405 0.3796441Ens indica que tenim punts de canvi als mesos de 1-2016, 8-2016, 1-2017, 8-2017, 1-2018, 8-2018, 1-2019, 8-2019. És a dir segueixen el mateix patró que abans.
Ara el valor de \(R^2\) de la segmentada és
summary(serie2_mall.seg)$adj.r.squared## [1] 0.8316885Aquest és més proper a 1 i podem assegurar que és bo.
Anem a visualitzar la regressió segmentada damunt les nostres dades
p_serie2_mall <- ggplot(mall[1:51,], aes(x=x, y=y)) + geom_line()+#dibuixam les nostres dades en línia contínua +
ggtitle('Regressió lineal i segmentada sobre les dades de Mallorca \ndurant la COVID') + xlab('índex del mes')+ ylab('Despeses mensuals en €')+ ylim(c(0,1300))
my.coef2_mall <- coef(serie2_mall.lm) #coeficients de la recta de regressió lineal
p_serie2_mall <- p_serie2_mall + geom_abline(intercept=my.coef2_mall[1], slope=my.coef2_mall[2], color="green") #afegim la recta de regressió lineal
my.model2_mall <- data.frame(x=1:51, y=fitted(serie2_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie2_mall <- p_serie2_mall + geom_line(data=my.model2_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines2_mall <- serie2_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie2_mall <- p_serie2_mall + geom_vline(xintercept=my.lines2_mall, linetype="dashed") #línies verticals en els punts de canvi
p_serie2_mall <- p_serie2_mall + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
ggplotly(p_serie2_mall)Per poder escriure la funció necessitam les pendents i els punts d’intersecció amb l’eix OY, que ens ho donen les següents funcions:
#PENDENTS
slope(serie2_mall.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -38.707 21.3120 -1.8162 -82.065 4.6519
## slope2 40.282 9.0058 4.4729 21.959 58.6040
## slope3 -61.599 15.0700 -4.0877 -92.259 -30.9400
## slope4 62.878 11.3910 5.5197 39.702 86.0540
## slope5 -75.010 15.0700 -4.9776 -105.670 -44.3510
## slope6 49.277 9.0058 5.4717 30.954 67.5990
## slope7 -72.938 11.3910 -6.4028 -96.114 -49.7610
## slope8 66.757 11.3910 5.8603 43.581 89.9340
## slope9 -85.307 15.0700 -5.6609 -115.970 -54.6480#PUNTS D'INTERSECCIÓ
intercept(serie2_mall.seg)## $x
## Est.
## intercept1 934.94
## intercept2 598.16
## intercept3 1718.80
## intercept4 -326.36
## intercept5 2831.40
## intercept6 -624.46
## intercept7 3641.70
## intercept8 -1978.00
## intercept9 5157.30Aleshores tenim
\(\hat{y}= \left\{ \begin{array}{lcc} -38.707x + 934.94 & si & x \leq \psi_1 \\ \\ 40.282x + 598.16 & si & \psi_1 < x \leq \psi_2 \\ \\ -61.599x + 1718.8 & si & \psi_2 < x \leq \psi_3 \\ \\ 62.878x - 326.36 & si & \psi_3 < x \leq \psi_4 \\ \\ -75.010x + 2831.4 & si & \psi_4 < x \leq \psi_5 \\ \\ 49.277x - 624.46 & si & \psi_5 < x \leq \psi_6 \\ \\ -72.938x + 3641.7 & si & \psi_6 < x \leq \psi_7\\ \\ 66.757x - 1978 & si & \psi_7 < x \leq \psi_8\\ \\ -85.307x + 5157.3 & si & \psi_8 < x \\ \end{array} \right.\)
on \(\psi_1= 4.26\), \(\psi_2= 11\), \(\psi_3= 16.43\), \(\psi_4= 22.9\), \(\psi_5= 27.8\), \(\psi_6= 34.91\), \(\psi_7= 40.23\) i \(\psi_8= 46.92\).
L’error del model de regressió segmentada és (ens interessa RMSE):
accuracy(serie2_mall.seg)## ME RMSE MAE MPE MAPE MASE
## Training set 4.458103e-15 38.33289 32.93155 -0.1657412 3.661674 0.3423069És a dir 38.33289 euros.
Anem a fer la predicció, en aquest cas també tenim que els següents punts de canvi es donaran en gener, agost i gener.
- La pendent de 8-2019 a 1-2020 i de 8-2020 a 1-2021 serà: (-61.599 - 75.010 -72.938)/3 = -69.849
- La pendent de 1-2020 a 8-2020 serà : (40.282 + 62.878 + 49.277 + 66.757)/4 = 54.799
Seguint el mateix procediment que abans obtenim que \(n_1 = 4431.55\), \(n_2 = -2050.14\) i \(n_3 = 5304.09\).
És a dir la funció per a la predicció de 2020 és :
\(\hat{y}= \left\{ \begin{array}{lcc} -69.849x + 4431.55 & si & x \leq \psi_9 \\ \\ 54.799x - 2050.14 & si & \psi_9 < x \leq \psi_{10} \\ \\ -69.849x + 5304.09 & si & \psi_{10} < x \\ \end{array} \right.\)
on \(\psi_9 = 52\) i \(\psi_8 = 59\).
#Ho graficam
p2_serie_mall <- ggplot(mall, aes(x=x, y=y)) + geom_line()+ #dibuixam les nostres dades en línia contínua
ggtitle('Predicció de Mallorca amb el model de regressió segmentada \ndurant la COVID') +
xlab('Índex del mes') +
ylab('Despeses mensuals en €') +
ylim(c(0,1300)) #dibuixam les nostres dades en línia contínua
my.model2_mall <- data.frame(x=1:51, y=fitted(serie2_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada
p2_serie_mall <- p2_serie_mall + geom_line(data=my.model2_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines2_mall <- serie2_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats
p2_serie_mall <- p2_serie_mall + geom_vline(xintercept=my.lines2_mall, linetype="dashed") #línies verticals en els punts de canvi
p2_serie_mall <- p2_serie_mall + geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
p2_serie_mall <- p2_serie_mall + geom_abline(intercept = 4431.55, slope=-69.849, colour="green") +
geom_abline(intercept = -2050.14 , slope=54.799, colour="blue") +
geom_abline(intercept = 5304.09 , slope=-69.849, colour="orange")
ggplotly(p2_serie_mall)Podem observar gràficament que per a l’any 2020 la predicció no és gens bona ja que hi va haver un fet inesperat que és la COVID.
Calculem aquest error:
o2_mall<-c(serie_mall[52:63]) # dades reals per fer predicció del 2020
v3=c(52:59)
f3 <- sapply(v3, function(x) 54.799*x-2050.14) #predicció de gener 2020 a agost 2020
v4=c(60:63)
f4 <- sapply(v4, function(x) -69.849*x + 5304.09) #predicció de setembre 2020 a desembre 2020
p2_mall <- c(f3,f4) #predicció de 2020
sqrt(sum((p2_mall-o2_mall)^2)/12)## [1] 442.1241En aquest cas l’error de la predicció amb aquest model és d’uns 425 euros.
- DECOMPOSE
Fem el mateix procediment que abans amb el mètode
decompose i predim l’any 2020.
La descomposició, en aquest cas, és la següent:
decompose_s2_mall <- decompose(serie2_mall)
plot(decompose_s2_mall,xlab="Any")
On les components de tendència són:
t2_mall <- decompose_s2_mall$trend #tendències de la sèrie sense el tros a predir abans del COVID
t2_mall## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 899.8458 899.0708 896.9637 895.5496 888.9975
## 2017 891.8125 896.2783 901.1758 905.1896 907.7217 908.4975 909.7354 914.6108
## 2018 924.3500 924.6162 923.1429 921.8554 921.3479 920.0171 918.4050 914.4658
## 2019 914.9921 918.6238 923.0717 925.9079 928.4362 930.7558 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 885.1633 886.3112 887.5750 888.9183
## 2017 920.9775 923.1212 922.9363 923.5758
## 2018 909.9675 910.4688 911.8138 912.6679
## 2019 NA NA NA NAI les estacionals:
s2_mall <- decompose_s2_mall$seasonal #estacionalitat
s2_mall <- s2_mall[4:15] #estacionalitat de gener a desembre
s2_mall## [1] -174.403157 -122.557740 -54.608435 -60.967983 3.885038 51.915663
## [7] 179.218371 198.616982 82.055593 71.194621 -116.073296 -58.275657Fem la predicció:
a2_mall <- c(s2_mall[7:12],s2_mall) #estacionalitat de juliol 2019 a desembre 2020 (es per poder fer la predicció)
pred2_decompose_mall <- sapply(a2_mall, function(x) 930.7558 + x) #predicció de juliol de 2019 a desembre 2020
p2_dec_mall<-c(serie2_mall[1:45], pred2_decompose_mall)
A2_mall<- data.frame("x" = mall[1:63,]$x, "y" = mall[1:63,]$y, "p"= p2_dec_mall)
grafica2_mall_dec <- ggplot(A2_mall)+
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+
labs(title="Predicció de Mallorca durant la COVID amb el model de descomposició", x="Índex del mes", y="Despesese mensuals en €")+
geom_vline(xintercept = 0)+ geom_hline(yintercept = 0)+ theme(panel.background = element_blank())
grafica2_mall_dec
Amb aquest mètode també s’observa que la caiguda del COVID no s’ha detectat, ja que de les dades on estam aprenent per fer la predicció no hi tenim valors pareguts als de la COVID.
L’error que es comet és:
sqrt(sum((c(serie_mall[52:63]- pred2_decompose_mall[7:18]))^2)/12)## [1] 390.5974Que també és un valor gros, tot i que no tant com amb el model anterior.
Vegem, igual que abans, la predicció amb la funció
predict():
prediccio2_mall <- predict(serie2_mall,h=12)
plot(prediccio2_mall, xlab="Any", ylab ="Despeses menusals en €")
Així pareix que ha de ser una bona predicció, però recordem que aquest cicle és el de la COVID, vegem com queda la predicció i la sèrie original:
df_prediccio2_mall <- data.frame(prediccio2_mall)
p2_ets_mall <- data.frame("x"= 1:63, "PointForecast"=c(serie2_mall,df_prediccio2_mall$Point.Forecast), "Lo80" = c(rep(NA,51),df_prediccio2_mall$Lo.80), "Hi80" = c(rep(NA,51),df_prediccio2_mall$Hi.80), "Lo95" = c(rep(NA,51),df_prediccio2_mall$Lo.95),"Hi95" = c(rep(NA,51),df_prediccio2_mall$Hi.95))
grafica_pred2_ets <- ggplot((mall[1:63,]))+
geom_ribbon(data = p2_ets_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p2_ets_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p2_ets_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(A,N,A) a Mallorca durant la COVID", y="Despeses mensuals en €", x="Índex del mes") +
geom_vline(xintercept = 0) + geom_hline(yintercept = 0)+ theme(panel.background = element_blank())
grafica_pred2_ets
Efectivament és una mala predicció, ja que no ha detectat la baixada de la COVID.
Si calculam l’error comés, aquest és de quasi 400 euros.
accuracy(prediccio2_mall,serie3_mall)## ME RMSE MAE MPE MAPE MASE
## Training set 1.768898 27.38848 22.10859 0.1014835 2.474142 0.6564259
## Test set -248.774469 392.70874 248.77447 -Inf Inf 7.3863591
## ACF1 Theil's U
## Training set 0.004436559 NA
## Test set 0.439341466 0- SARIMA
Quin model proposam nosaltres? Vegem l’ACF i el PACF.
par(mfrow=c(1,2))
acf(serie2_mall)
pacf(serie2_mall)
Per a la part regular, de l’ACF obtenim un MA(2) i del PACF un AR(1). Notem que no fa falta fer cap diferenciació a la part regular ja que no considram que hi hagi tendència. No obstant, sí que observam que hi ha una certa estacionalitat per la forma que tenen les barres en els correlogrames.
Llavors fem una diferenciació d’ordre 12:
serie2_mall_diff <- diff(serie2_mall,12)
plot(serie2_mall_diff, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")
Ara ja no s’observa l’estacionalitat, llavors hem de fer una diferenciació D=1. Vegem, ara, quins són els correlogrames
par(mfrow=c(1,2))
acf(serie2_mall_diff,lag=36)
pacf(serie2_mall_diff,lag=36)
Amb aquests obtenim un MA(1) i un AR(1) per a la part estacional. És a dir, per predir un més hem de tenir en compte el valor i els errors comesos d’aquest mateix mes un cicle anterior.
Llavors el model que nosaltres proposam es un ARIMA(1,0,2)(1,1,1)[12]
R proposa el següent model:
auto.arima(serie2_mall)## Series: serie2_mall
## ARIMA(0,0,0)(1,1,0)[12] with drift
##
## Coefficients:
## sar1 drift
## -0.7165 0.9022
## s.e. 0.1136 0.2674
##
## sigma^2 = 927: log likelihood = -191.85
## AIC=389.7 AICc=390.39 BIC=394.69Fixem-nos que és el mateix model que proposa R per a la serie1, és a dir, la que va d’octubre de 2015 a desembre de 2018.
Els dos models que tenim són:
model3_mall <- arima(serie2_mall, order=c(1,0,2), seasonal = list(order=c(1,1,1), period=12))
model4_mall <- arima(serie2_mall, order=c(0,0,0), seasonal = list( order=c(1,1,0), period=12))I els seus errors
accuracy(model3_mall)## ME RMSE MAE MPE MAPE MASE
## Training set 5.895861 22.21249 15.13338 0.5278811 1.683231 0.2080521
## ACF1
## Training set -0.02958489accuracy(model4_mall)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 10.85733 30.49463 21.65555 0.9581157 2.367599 0.2977183 0.2066374Les prediccions són:
fc_m3_mall <- forecast(model3_mall, h=12)
fc_m4_mall <- forecast(model4_mall,h=12)#PRIMER MODEL
pre2_mall <- data.frame("Point Forecast" = serie2_mall, "Lo 80" = rep(NA,51), "Hi 80"= rep(NA,51), "Lo 95" = rep(NA,51), "Hi 95" = rep(NA,51))
pred_m3_mall <-data.frame(fc_m3_mall)
grafica_m3_mall <- data.frame("x" = 1:63, "PointForecast" = c(pre2_mall$Point.Forecast,pred_m3_mall$Point.Forecast), "Lo80" = c(pre2_mall$Lo.80, pred_m3_mall$Lo.80), "Hi80"= c(pre2_mall$Hi.80, pred_m3_mall$Hi.80), "Lo95" = c(pre2_mall$Lo.95, pred_m3_mall$Lo.95), "Hi95" = c(pre2_mall$Hi.95, pred_m3_mall$Hi.95))
grafica3_mall <- ggplot(mall[1:63,])+
geom_ribbon(data = grafica_m3_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m3_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m3_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red") +
labs(title="Prediccó de Mallorca amb el model ARIMA (1,0,2)(1,1,1)[12] durant la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica3_mall
#SEGON MODEL
pred_m4_mall <-data.frame(fc_m4_mall)
grafica_m4_mall <- data.frame("x" = 1:63, "PointForecast" = c(pre2_mall$Point.Forecast,pred_m4_mall$Point.Forecast), "Lo80" = c(pre2_mall$Lo.80, pred_m4_mall$Lo.80), "Hi80"= c(pre2_mall$Hi.80, pred_m4_mall$Hi.80), "Lo95" = c(pre2_mall$Lo.95, pred_m4_mall$Lo.95), "Hi95" = c(pre2_mall$Hi.95, pred_m4_mall$Hi.95))
grafica4_mall <- ggplot(mall[1:63,])+
geom_ribbon(data = grafica_m4_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m4_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m4_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red") +
labs(title="Predicció de Mallorca amb el model ARIMA (0,0,0)(1,1,0)[12] durant la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica4_mall
Amb un error de 384.5 i 383.6:
accuracy(fc_m3_mall,serie_mall[52:63], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set 5.895861 22.21249 15.13338 0.5278811 1.683231 0.2080521
## Test set -243.110836 384.52595 243.11084 -Inf Inf 3.3422626
## ACF1
## Training set -0.02958489
## Test set NAaccuracy(fc_m4_mall,serie_mall[52:63], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set 10.85733 30.49463 21.65555 0.9581157 2.367599 0.2977183
## Test set -237.63235 383.58017 238.11170 -Inf Inf 3.2735350
## ACF1
## Training set 0.2066374
## Test set NAEn els dos models veiem que la predicció és dolenta ja que segueixen el mateix patró que els cicles anteriors i no preveuen la COVID. També podem observar que l’error que es comet és bastant alt.
Igual que abans, esdudiem amb més detall aquests errors:
checkresiduals(fc_m3_mall)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2)(1,1,1)[12]
## Q* = 4.7827, df = 5, p-value = 0.443
##
## Model df: 5. Total lags used: 10checkresiduals(fc_m4_mall)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,0)(1,1,0)[12]
## Q* = 13.035, df = 9, p-value = 0.161
##
## Model df: 1. Total lags used: 10par(mfrow=c(1,2))
qqPlot(fc_m3_mall$residuals, id=FALSE, mean=mean(fc_m3_mall$residuals), sd=sd(fc_m3_mall$residuals))
qqPlot(fc_m4_mall$residuals, id=FALSE, mean=mean(fc_m4_mall$residuals), sd=sd(fc_m4_mall$residuals))
shapiro.test(fc_m3_mall$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m3_mall$residuals
## W = 0.93903, p-value = 0.01123shapiro.test(fc_m4_mall$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m4_mall$residuals
## W = 0.93172, p-value = 0.005797lillie.test(fc_m3_mall$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m3_mall$residuals
## D = 0.14786, p-value = 0.00709lillie.test(fc_m4_mall$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m4_mall$residuals
## D = 0.20836, p-value = 7.9e-06Obtenim que els tests de normalitat ens donen uns p-valors molt petits. No obstant, en els dos models els residus no segueixen cap estructura i les barres dels ACF es mantenen dins les bandes de confiança. També observam que, als qqPlot, tot i haver-hi valors fora de la regió blava i de no estar sobre la diagonal, la majoria de valors es troben al voltant de 0.
Vegem un resum dels errors que cometen cada un dels models en aquest cas:
| reg. segmentada | descomposició clàssica | ETS(A,N,A) | ARIMA (1,0,2)(1,1,1)[12] | ARIMA (0,0,0)(1,1,0)[12] | |
|---|---|---|---|---|---|
| error model | 38.33 | 27.39 | 22.21 | 30.49 | |
| error predicció | 424.78 | 390.6 | 392.71 | 384.53 | 383.58 |
2.3. DESPRÉS DE LA COVID
- REGRESSIÓ SEGMENTADA
El valor de \(R^2\) per a la regressió lineal és
serie3_mall.lm <- lm(y~x, data= mall[1:63,])
summary(serie3_mall.lm)$adj.r.squared## [1] 0.02681106Vegem els punts de canvi:
punts_canvi_serie3_mall <-selgmented(serie3_mall.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 ..
## (search truncated at 6 breakpoints due to increasing values of BIC)
##
## BIC to detect no. of breakpoints:
## 0 1 2 3 4 5 6
## 670.4792 673.4233 681.6172 688.4189 695.2176 701.0818 706.7486
##
## No. of selected breakpoints: 0Fixem-nos que si demanam a R que ens indiqui quants de punts de canvi estima amb la sèrie completa, és a dir, amb la COVID, no detecta cap punt de canvi. Aquest fet s’eslica amb detall al treball.
serie3_mall.seg <- segmented(serie3_mall.lm, seg.Z=~x, psi = c(5,10,17,23,28,35,40,47,55,58))
summary(serie3_mall.seg)$psi## Initial Est. St.Err
## psi1.x 5 5.783722 1.7320089
## psi2.x 10 10.400505 1.1467103
## psi3.x 17 16.853561 1.0053175
## psi4.x 23 22.744051 0.8995765
## psi5.x 28 27.944118 1.0394421
## psi6.x 35 34.740660 0.9251495
## psi7.x 40 40.250732 0.7944349
## psi8.x 47 47.223848 0.5600066
## psi9.x 55 56.000020 0.4535518
## psi10.x 58 58.017477 0.3070058Aquests es donen als mesos següents: 3-2026, 7-2016, 2-2017, 8-2017, 1-2018, 8-2018, 1-2019, 8-2019, 5-2020, 7-2020.
Ara, el valor de l’\(R^2\) de la segmentada és
summary(serie3_mall.seg)$adj.r.squared## [1] 0.7222148Anem a visualitzar la regressió segmentada damunt les nostres dades
p_serie3_mall <- ggplot(mall[1:63,], aes(x=x, y=y)) + geom_line()+ #dibuixam les nostres dades en línia contínua
ggtitle('Regressió lineal i segmentada sobre les dades \nde Mallorca després de la COVID') + xlab('índex del mes')+ ylab('Despeses mensuals en €')+
ylim(c(0,1300))
#dibuixam les nostres dades en línia contínua
my.coef3_mall <- coef(serie3_mall.lm) #coeficients de la recta de regressió lineal
p_serie3_mall <- p_serie3_mall + geom_abline(intercept=my.coef3_mall[1], slope=my.coef3_mall[2], color="green") #afegim la recta de regressió lineal
my.model3_mall <- data.frame(x=1:63, y=fitted(serie3_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie3_mall <- p_serie3_mall + geom_line(data=my.model3_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines3_mall <- serie3_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie3_mall <- p_serie3_mall + geom_vline(xintercept=my.lines3_mall, linetype="dashed") #línies verticals en els punts de canvi
p_serie3_mall <- p_serie3_mall + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
ggplotly(p_serie3_mall)Per escriure la funció a trossos tenim:
#PUNTS D'INTERSECCIÓ
intercept(serie3_mall.seg)## $x
## Est.
## intercept1 897.69
## intercept2 449.26
## intercept3 1601.80
## intercept4 -444.27
## intercept5 2752.50
## intercept6 -722.00
## intercept7 3636.00
## intercept8 -2243.00
## intercept9 6553.70
## intercept10 -23235.00
## intercept11 5426.40#PENDENTS
slope(serie3_mall.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -19.420 32.900 -0.59027 -85.863 47.0230
## slope2 58.115 32.900 1.76640 -8.328 124.5600
## slope3 -52.703 24.870 -2.11910 -102.930 -2.4764
## slope4 68.701 24.870 2.76240 18.475 118.9300
## slope5 -71.853 32.900 -2.18400 -138.300 -5.4096
## slope6 52.484 19.662 2.66940 12.777 92.1920
## slope7 -72.961 24.870 -2.93370 -123.190 -22.7340
## slope8 73.101 19.662 3.71800 33.394 112.8100
## slope9 -113.180 13.431 -8.42630 -140.300 -86.0520
## slope10 418.770 147.130 2.84620 121.630 715.9100
## slope11 -75.248 32.900 -2.28720 -141.690 -8.8050L’error de la regressió és:
accuracy(serie3_mall.seg)## ME RMSE MAE MPE MAPE MASE
## Training set 3.609932e-15 83.93007 53.6569 -Inf Inf 0.4360175Anem a fer la predicció per a l’any 2021. Recordem que els punts de canvi es donen el 3-2026, 7-2016, 2-2017, 8-2017, 1-2018, 8-2018, 1-2019, 8-2019, 5-2020, 7-2020. Degut a la pertorbació del COVID, sí que hi segueix havent un punt de canvi en estiu i l’altre en hivern successivament però ara no és donen al mateix mes. Per això, per predir els següents punts de canvi agafarem el mes més freqüent. Aleshores els següents punts de canvi seran en 1-2021, 8-2021 i 1-2022.
- La pendent de 7-2020 a 1-2021 i de 8-2021 a 1-2022 serà : -77.674
- La pendent de 1-2021 a 8-2021 serà : 134.234
Els nous punts d’intersecció es calculen de forma anàloga que als casos anteriors. Per a les tres noves rectes obtenim que són \(n_1=5563.7164\), \(n_2=-7998.396\) i \(n_3=7047.07\).
Aleshores, la predicció de 2021 serà:
\(\hat{y}= \left\{ \begin{array}{lcc} -77.674x + 5563.716 & si & x \leq \psi_{11} \\ \\ 134.234x - 7998.396 & si & \psi_{11} < x \leq \psi_{12} \\ \\ -77.674x + 7047.07 & si & \psi_{12} < x \\ \end{array} \right.\)
on \(\psi_7 = 64\) i \(\psi_8 = 71\).
#visualitzem-la
p3_serie_mall <- ggplot(mall, aes(x=x, y=y)) + geom_line()+ #dibuixam les nostres dades en línia contínua
ggtitle('Predicció de Mallorca amb el model de \nregressió segmentada després de la COVID') + xlab('índex del mes')+ ylab('Despeses mensuals en €')+
ylim(c(0,1600))
my.model3_mall <- data.frame(x=1:63, y=fitted(serie3_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada
p3_serie_mall <- p3_serie_mall + geom_line(data=my.model3_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines3_mall <- serie3_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats
p3_serie_mall <- p3_serie_mall + geom_vline(xintercept=my.lines3_mall, linetype="dashed") #línies verticals en els punts de canvi
p3_serie_mall <- p3_serie_mall + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0)+theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
p3_serie_mall <- p3_serie_mall + geom_abline(intercept = 5563.7164, slope=-77.674, colour="green") +
geom_abline(intercept = -7998.396 , slope= 134.234, colour="blue") +
geom_abline(intercept = 7047.07 , slope=-77.674, colour="orange")
ggplotly(p3_serie_mall)A simple vista observam que la predicció no és massa bona, ja que la predicció que es prediu pel mes d’estiu de l’any 2021 és d’uns 1500€ aproximadament, mentre que el valor real és d’uns 1100€.
Vegem quin és aquest error que es comet:
o3_mall<-c(serie_mall[64:75]) # dades reals per fer predicció del 2021
v5=c(64:71)
f5 <- sapply(v5, function(x) 134.234*x-7998.396) #predicció de gener 2021 a agost 2021
v6=c(72:75)
f6 <- sapply(v6, function(x) -77.674*x + 7047.07) #predicció de setembre 2021 a desembre 2021
p3_mall <- c(f5,f6) #predicció de 2021
sqrt(sum((p3_mall-o3_mall)^2)/12)## [1] 269.9865- DECOMPOSE
Visualitzem la sèrie descomposada:
decompose_s3_mall <- decompose(serie3_mall)
plot(decompose_s3_mall, xlab="Any")
Els components de tendència són:
t3_mall <- decompose_s3_mall$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t3_mall## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 899.8458 899.0708 896.9637 895.5496 888.9975
## 2017 891.8125 896.2783 901.1758 905.1896 907.7217 908.4975 909.7354 914.6108
## 2018 924.3500 924.6162 923.1429 921.8554 921.3479 920.0171 918.4050 914.4658
## 2019 914.9921 918.6238 923.0717 925.9079 928.4362 930.7558 931.6700 934.0650
## 2020 750.7167 737.9196 727.9308 713.3075 700.2629 692.3146 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 885.1633 886.3112 887.5750 888.9183
## 2017 920.9775 923.1212 922.9363 923.5758
## 2018 909.9675 910.4688 911.8138 912.6679
## 2019 931.8875 891.1542 816.3763 767.6412
## 2020 NA NA NA NAI els d’estacionalitat:
s3_mall <- decompose_s3_mall$seasonal #estacionalitat
s3_mall <- s3_mall[4:15] #estacionalitat de gener a desembre
s3_mall## [1] -127.84021 -68.39937 -27.68521 -186.60614 -132.11481 55.96136
## [7] 191.44261 213.85281 92.39615 86.61625 -78.02771 -19.59573Anem a fer la predicció:
a3_mall <- c(s3_mall[7:12],s3_mall) #estacionalitat de juliol 2020 a desembre 2021 (es per poder fer la predicció)
pred3_decompose_mall <- sapply(a3_mall, function(x) 692.3146 + x) #predicció de juliol de 2019 a desembre 2020
p3_dec_mall<-c(serie3_mall[1:57], pred3_decompose_mall)
A3_mall<- data.frame("x" = mall[1:75,]$x, "y" = mall[1:75,]$y, "p"= p3_dec_mall)
grafica3_mall_dec <- ggplot(A3_mall)+
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+
labs(title="Predicció de Mallorca després de la COVID amb el model de descomposició", x="Índex del mes", y="Despeses mensualns en €")+geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica3_mall_dec
Observam que la predicció, en vermell, segueix més o menys el mateix patró que la sèrie original. No obstant això, es troba a un nivell més baix degut a les dades anteriors que té per aprendre, on hi ha l’efecte de la COVID.
Calculem l’error de la predicció:
sqrt(sum((c(serie_mall[64:75]- pred3_decompose_mall[7:18]))^2)/12)## [1] 280.1672Amb la funció predict() la predicció seria la
següent:
prediccio3_mall <- predict(serie3_mall,h=12)
plot(prediccio3_mall, xlab="Any", ylab="Despeses mensuals en €")
Observam que, a diferència de les series anteriors amb aquesta funció, sense la COVID, la predicció és molt dolenta. Vegem-la juntament amb les nostres dades:
df_prediccio3_mall <- data.frame(prediccio3_mall)
p3_ets_mall <- data.frame("x"= 1:75, "PointForecast"=c(serie3_mall,df_prediccio3_mall$Point.Forecast), "Lo80" = c(rep(NA,63),df_prediccio3_mall$Lo.80), "Hi80" = c(rep(NA,63),df_prediccio3_mall$Hi.80), "Lo95" = c(rep(NA,63),df_prediccio3_mall$Lo.95),"Hi95" = c(rep(NA,63),df_prediccio3_mall$Hi.95))
grafica_pred3_ets <- ggplot((mall[1:75,]))+
geom_ribbon(data = p3_ets_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p3_ets_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p3_ets_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(A,N,N) a Mallorca després de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica_pred3_ets
Aquest fet és degut a que el model ETS dóna més pes a les observacions anteriors properes i no a les llunyanes, i no està considerant estacionalitat degut a la detecció de la COVID.
accuracy(prediccio3_mall, serie_mall)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -3.019093 162.3973 95.90263 -Inf Inf 1.119973 0.1219558
## Test set 251.816426 282.9966 251.81643 24.85994 24.85994 2.940770 0.6915062
## Theil's U
## Training set NA
## Test set 3.250169- SARIMA
Quin model proposam nosaltres? Vegem l’ACF i el PACF:
par(mfrow=c(1,2))
acf(serie3_mall, lag.max=36)
pacf(serie3_mall, lag.max=36)
Per a la part regular obtenim un MA(1), un AR(2). Observam una clara estacionalitat el l’ACF però com les barres no surten no la consideram significativa. Llavors d=1.
Diferenciem la part estacional i calculem el seu ACF i PACF:
serie3_mall_diff <- diff(serie3_mall,12)
plot.ts(serie3_mall_diff, main="Sèrie sense estacionalitat", ylab="Sèrie diferenciada", xlab="Any")
par(mfrow=c(1,2))
acf(serie3_mall_diff)
pacf(serie3_mall_diff)
Fent una diferenciació d’ordre 12, per llevar estacionalitat, i tornant a calcular els ACF i PACF obtenim: que l’ordre del model autoregressiu és P=2 i el de mitjana mòbil Q=1.
Aleshores la nostra proposta és un model ARIMA(2,1,1)(2,0,1)[12]
R proposa el següent model:
auto.arima(serie3_mall)## Series: serie3_mall
## ARIMA(1,0,1)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ma1 sar1 mean
## 0.3171 0.7225 0.6045 831.5038
## s.e. 0.1369 0.0994 0.1267 75.8692
##
## sigma^2 = 14934: log likelihood = -393.43
## AIC=796.86 AICc=797.91 BIC=807.58Llavors tenim els següents models:
model5_mall <- arima(serie3_mall, order=c(2,1,1), seasonal = list( order=c(2,0,1), period=12))
model6_mall <- arima(serie3_mall, order=c(1,0,1), seasonal = list( order=c(1,0,0), period=12))Amb un error de 112.18 i 118.26, respectivament:
accuracy(model5_mall)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -3.612028 127.6573 76.78002 NaN Inf 0.7880336 -0.03300271accuracy(model6_mall)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -4.030176 118.2628 71.00525 -Inf Inf 0.7287641 0.0005558812Les prediccions d’aquests models són:
fc_m5_mall <- forecast(model5_mall, h=12)
fc_m6_mall <- forecast(model6_mall,h=12)#grafiquem-les
#primer model
pre3_mall <- data.frame("Point Forecast" = serie3_mall, "Lo 80" = rep(NA,63), "Hi 80"= rep(NA,63), "Lo 95" = rep(NA,63), "Hi 95" = rep(NA,63))
pred_m5_mall <-data.frame(fc_m5_mall)
grafica_m5_mall <- data.frame("x" = 1:75, "PointForecast" = c(pre3_mall$Point.Forecast,pred_m5_mall$Point.Forecast), "Lo80" = c(pre3_mall$Lo.80, pred_m5_mall$Lo.80), "Hi80"= c(pre3_mall$Hi.80, pred_m5_mall$Hi.80), "Lo95" = c(pre3_mall$Lo.95, pred_m5_mall$Lo.95), "Hi95" = c(pre3_mall$Hi.95, pred_m5_mall$Hi.95))
grafica5_mall <- ggplot(mall[1:75,])+
geom_ribbon(data = grafica_m5_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m5_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m5_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció de Mallorca amb el model ARIMA (2,1,1)(2,0,1)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica5_mall
#segon model
pred_m6_mall <-data.frame(fc_m6_mall)
grafica_m6_mall <- data.frame("x" = 1:75, "PointForecast" = c(pre3_mall$Point.Forecast,pred_m6_mall$Point.Forecast), "Lo80" = c(pre3_mall$Lo.80, pred_m6_mall$Lo.80), "Hi80"= c(pre3_mall$Hi.80, pred_m6_mall$Hi.80), "Lo95" = c(pre3_mall$Lo.95, pred_m6_mall$Lo.95), "Hi95" = c(pre3_mall$Hi.95, pred_m6_mall$Hi.95))
grafica6_mall <- ggplot(mall[1:75,])+
geom_ribbon(data = grafica_m6_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m6_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m6_mall, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció de Mallorca amb el model ARIMA (1,0,1)(1,0,0)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica6_mall
Tant a d’un model com a l’altre observam que la predicció no és molt bona. De la mateixa manera que als models anteriors, aquest fet és degut a l’alteració de les dades per la COVID.
Vegem quin és l’error de la predicció per a cada model
accuracy(fc_m5_mall,serie_mall[64:75], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set -3.612028 127.6573 76.78002 NaN Inf 0.7880336
## Test set 422.131202 475.4906 422.13120 43.52191 43.52191 4.3325538
## ACF1
## Training set -0.03300271
## Test set NAaccuracy(fc_m6_mall,serie_mall[64:75], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set -4.030176 118.2628 71.00525 -Inf Inf 0.7287641
## Test set 236.006237 308.5205 240.53442 24.01452 24.63381 2.4687309
## ACF1
## Training set 0.0005558812
## Test set NAObtenim que l’error és de 475.5 i 308.5.
Estudiem-los:
checkresiduals(fc_m5_mall)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,1)(2,0,1)[12]
## Q* = 9.1059, df = 7, p-value = 0.2451
##
## Model df: 6. Total lags used: 13checkresiduals(fc_m6_mall)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1)(1,0,0)[12] with non-zero mean
## Q* = 4.3222, df = 10, p-value = 0.9316
##
## Model df: 3. Total lags used: 13par(mfrow=c(1,2))
qqPlot(fc_m5_mall$residuals, id=FALSE, mean=mean(fc_m5_mall$residuals), sd=sd(fc_m5_mall$residuals))
qqPlot(fc_m6_mall$residuals, id=FALSE, mean=mean(fc_m6_mall$residuals), sd=sd(fc_m6_mall$residuals))
shapiro.test(fc_m5_mall$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m5_mall$residuals
## W = 0.79747, p-value = 6.838e-08shapiro.test(fc_m6_mall$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m6_mall$residuals
## W = 0.69775, p-value = 4.259e-10lillie.test(fc_m5_mall$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m5_mall$residuals
## D = 0.18624, p-value = 1.068e-05lillie.test(fc_m6_mall$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m6_mall$residuals
## D = 0.16427, p-value = 0.0002185De la mateixa manera que als casos anteriors, els errors no segueixen una normal.
Vegem un resum dels errors de tots els models per aquest cas:
| reg. segmentada | descomposició clàssica | ETS(A,N,N) | ARIMA (2,1,1)(2,0,1)[12] | ARIMA (1,0,1)(1,0,0)[12] | |
|---|---|---|---|---|---|
| error model | 83.93 | 162.4 | 112.18 | 118.26 | |
| error predicció | 269.987 | 280.17 | 283 | 451.8 | 308.52 |